MODULE fefmm_LIP3_Fij

   USE fmm_global_paras
   IMPLICIT NONE
   PRIVATE
   ! Public procedures
   PUBLIC :: fefmm_build_Fij_LIP3

   REAL(REALK),   PARAMETER :: T_VALUE_THR = 0.5D0
   INTEGER(INTK), PARAMETER :: LIPN = 3

CONTAINS

!-------------------------------------------------------------------------------

   SUBROUTINE fefmm_build_Fij_LIP3(FEdim,NFS,FEgrain,t,Fij)

      IMPLICIT NONE
      INTEGER(INTK), INTENT(IN)  :: FEdim, NFS
      REAL(REALK),   INTENT(IN)  :: FEgrain, t
      REAL(REALK),   INTENT(OUT) :: Fij(:,:)

      INTEGER(INTK) :: i,k, u,v

      Fij = zero

      DO i = 1, FEdim 
         v = i + FEdim*((NFS-1)/2)  ! V-point

         DO k = 1, NFS*( (FEdim -1)/LIPN )   ! finite elements
            u = 1 + LIPN*(k-1)    ! basic rho-point

            ! 1st LIP
            Fij(u,i)   = Fij(u,i)   + get_F1_LIP3(t,FEgrain,v,u)
            ! 2nd LIP
            Fij(u+1,i) = Fij(u+1,i) + get_F2_LIP3(t,FEgrain,v,u)
            ! 3rd LIP
            Fij(u+2,i) = Fij(u+2,i) + get_F3_LIP3(t,FEgrain,v,u)
            ! 4th LIP
            Fij(u+3,i) = Fij(u+3,i) + get_F4_LIP3(t,FEgrain,v,u)
!print *, u,i, v-u, Fij(u,i)

         END DO

      END DO

   END SUBROUTINE fefmm_build_Fij_LIP3

!-------------------------------------------------------------------------------

   FUNCTION get_F1_LIP3(t,lambda,v,u)

      IMPLICIT NONE
      REAL(REALK),   INTENT(IN) :: t
      REAL(REALK),   INTENT(IN) :: lambda
      INTEGER(INTK), INTENT(IN) :: v  ! this is the V-point
      INTEGER(INTK), INTENT(IN) :: u  ! this is the rho-point
      REAL(REALK) :: get_F1_LIP3

      INTEGER(INTK) :: k
      REAL(REALK)   :: F1, F2

      k = v-u

      IF (ABS(t) < T_VALUE_THR) THEN

      get_F1_LIP3 = (0.21D2 / 0.80D2 * lambda ** 9 * k ** 6 - 0.27D2    &
        / 0.20D2 * lambda ** 9 * k ** 5 - lambda ** 9 * k ** 7 / 0.20D2    &
        + k ** 8 * lambda ** 9 / 0.64D2 - 0.4617D4 / 0.440D3 * lambda ** 9    &
        * k + 0.2187D4 / 0.160D3 * lambda ** 9 * k ** 2 - 0.81D2 / 0.8D1    &
        * lambda ** 9 * k ** 3 + 0.297D3 / 0.64D2 * lambda ** 9 * k ** 4    &
        + 0.12393D5 / 0.3520D4 * lambda ** 9) * t ** 8 + (-0.9D1 / 0.16D2    &
        * lambda ** 7 * k ** 4 + 0.3D1 / 0.20D2 * lambda ** 7 * k ** 5    &
        - k ** 6 * lambda ** 7 / 0.16D2 - 0.2187D4 / 0.1120D4    &
        * lambda ** 7 - 0.891D3 / 0.224D3 * lambda ** 7 * k ** 2    &
        + 0.27D2 / 0.14D2 * lambda ** 7 * k ** 3 + 0.243D3 / 0.56D2    &
        * lambda ** 7 * k) * t ** 6 + (-0.81D2 / 0.70D2 * lambda ** 5    &
        * k + 0.891D3 / 0.1120D4 * lambda ** 5 + 0.27D2 / 0.40D2    &
        * lambda ** 5 * k ** 2 - 0.3D1 / 0.10D2 * lambda ** 5 * k ** 3    &
        + 0.3D1 / 0.16D2 * k ** 4 * lambda ** 5) * t ** 4 + (-0.9D1    &
        / 0.40D2 * lambda ** 3 - 0.3D1 / 0.8D1 * k ** 2 * lambda ** 3    &
        + 0.3D1 / 0.10D2 * lambda ** 3 * k) * t ** 2 + 0.3D1    &
        / 0.8D1 * lambda

      ELSE

      F1 = (-0.11D2 * lambda ** 2 * t ** 2 * exp(-t ** 2 * k ** 2         & 
         * lambda ** 2) - t ** 2 * k ** 2 * lambda ** 2 * exp(-t ** 2 *  & 
         k ** 2 * lambda ** 2) + 0.6D1 * lambda ** 2 * t ** 2 * k    &
         * exp(-t ** 2 * k ** 2 * lambda ** 2) - exp(-t ** 2 * k ** 2    &
         * lambda ** 2) - 0.3D1 * exp(-lambda ** 2 * t ** 2    &
         * (k - 0.3D1) ** 2) * k * t ** 2 * lambda ** 2 + 0.2D1    &
         * exp(-lambda ** 2 * t ** 2 * (k - 0.3D1) ** 2) * t ** 2    &
         * lambda ** 2 + exp(-lambda ** 2 * t ** 2 * (k - 0.3D1) ** 2)    &
         * k ** 2 * t ** 2 * lambda ** 2 + exp(-lambda ** 2 * t ** 2    &
         * (k - 0.3D1) ** 2)) / t ** 4 / lambda ** 3 / 0.12D2

      F2 = (((k ** 2 / 0.2D1 - k ** 3 / 0.12D2 - 0.11D2 / 0.12D2 * k    &
         + 0.1D1 / 0.2D1) * erf(t * k * lambda) + k ** 3 * erf(-0.3D1    &
         * lambda * t + t * k * lambda) / 0.12D2 + 0.11D2 / 0.12D2 * k    &
         * erf(-0.3D1 * lambda * t + t * k * lambda) - erf(-0.3D1 * lambda    &
         * t + t * k * lambda) / 0.2D1 - k ** 2 * erf(-0.3D1 * lambda * t    &
         + t * k * lambda) / 0.2D1) / t + ((-0.1D1 / lambda ** 2 * k    &
         / 0.8D1 + 0.1D1 / lambda ** 2 / 0.4D1) * erf(t * k * lambda)    &
         + 0.1D1 / lambda ** 2 * k * erf(-0.3D1 * lambda * t + t * k    &
         * lambda) / 0.8D1 - 0.1D1 / lambda ** 2 * erf(-0.3D1 * lambda    &
         * t + t * k * lambda) / 0.4D1) / t ** 3) * sqrt(0.3141592654D1)

      get_F1_LIP3 = F1 + F2

      END IF

   END FUNCTION get_F1_LIP3

!-------------------------------------------------------------------------------

   FUNCTION get_F2_LIP3(t,lambda,v,u)

      IMPLICIT NONE
      REAL(REALK),   INTENT(IN) :: t
      REAL(REALK),   INTENT(IN) :: lambda
      INTEGER(INTK), INTENT(IN) :: v  ! this is the V-point
      INTEGER(INTK), INTENT(IN) :: u  ! this is the rho-point
      REAL(REALK) :: get_F2_LIP3

      INTEGER(INTK) :: k
      REAL(REALK)   :: F1, F2

      k = v-u

      IF (ABS(t) < T_VALUE_THR) THEN

      get_F2_LIP3 = (-0.5103D4 / 0.80D2 * lambda ** 9 * k ** 2   &
        + 0.729D3 / 0.16D2 * lambda ** 9 * k ** 3 - 0.1215D4 / 0.64D2   &
        * lambda ** 9 * k ** 4 + 0.3D1 / 0.64D2 * lambda ** 9 * k ** 8   &
        - 0.9D1 / 0.40D2 * lambda ** 9 * k ** 7 + 0.81D2 / 0.20D2   &
        * lambda ** 9 * k ** 5 + 0.2187D4 / 0.44D2 * lambda ** 9   &
        * k - 0.59049D5 / 0.3520D4 * lambda ** 9) * t ** 8   &
        + (-0.3D1 / 0.16D2 * lambda ** 7 * k ** 6 + 0.729D3 / 0.80D2   &
        * lambda ** 7 - 0.81D2 / 0.14D2 * lambda ** 7 * k ** 3 + 0.3645D4   &
        / 0.224D3 * lambda ** 7 * k ** 2 - 0.2187D4 / 0.112D3   &
        * lambda ** 7 * k + 0.27D2 / 0.40D2 * lambda ** 7 * k ** 5)   &
        * t ** 6 + (0.243D3 / 0.70D2 * lambda ** 5 * k - 0.729D3   &
        / 0.224D3 * lambda ** 5 - 0.27D2 / 0.20D2 * lambda ** 5 * k ** 3   &
        + 0.9D1 / 0.16D2 * lambda ** 5 * k ** 4) * t ** 4 + (0.27D2   &
        / 0.20D2 * lambda ** 3 * k - 0.9D1 / 0.8D1 * lambda ** 3 * k ** 2)   &
        * t ** 2 + 0.9D1 / 0.8D1 * lambda

      ELSE

      F1 = -(-0.2D1 * exp(-dble(lambda ** 2 * t ** 2 * (k - 3) ** 2))   &
         * dble(k) * dble(t ** 2) * dble(lambda ** 2) - 0.6D1   &
         * dble(lambda ** 2) * dble(t ** 2) * exp(-dble(t ** 2 * k ** 2   &
         * lambda ** 2)) - dble(t ** 2) * dble(k ** 2) * dble(lambda ** 2)   &
         * exp(-dble(t ** 2 * k ** 2 * lambda ** 2)) + 0.5D1   &
         * dble(lambda ** 2) * dble(t ** 2) * dble(k) * exp(-dble(t ** 2   &
         * k ** 2 * lambda ** 2)) + exp(-dble(lambda ** 2 * t ** 2   &
         * (k - 3) ** 2)) + exp(-dble(lambda ** 2 * t ** 2 * (k - 3) ** 2))   &
         * dble(k ** 2) * dble(t ** 2) * dble(lambda ** 2)   &
         - exp(-dble(t ** 2 * k ** 2 * lambda ** 2))) / dble(t ** 4)   &
         / dble(lambda ** 3) / 0.4D1

      F2 = (((-0.5D1 / 0.4D1 * k ** 2 + k ** 3 / 0.4D1 + 0.3D1   &
         / 0.2D1 * k) * erf(t * k * lambda) + 0.5D1 / 0.4D1 * k ** 2   &
         * erf(-0.3D1 * lambda * t + t * k * lambda) - k ** 3 * erf(-0.3D1   &
         * lambda * t + t * k * lambda) / 0.4D1 - 0.3D1 / 0.2D1 * k   &
         * erf(-0.3D1 * lambda * t + t * k * lambda)) / t + ((-0.5D1   &
         / 0.8D1 / lambda ** 2 + 0.3D1 / 0.8D1 / lambda ** 2 * k)   &
         * erf(t * k * lambda) + 0.5D1 / 0.8D1 / lambda ** 2 * erf(-0.3D1   &
         * lambda * t + t * k * lambda) - 0.3D1 / 0.8D1 / lambda ** 2 * k   &
         * erf(-0.3D1 * lambda * t + t * k * lambda)) / t ** 3)   &
         * sqrt(0.3141592654D1)

      get_F2_LIP3 = F1 + F2

      END IF

   END FUNCTION get_F2_LIP3

!-------------------------------------------------------------------------------

   FUNCTION get_F3_LIP3(t,lambda,v,u)

      IMPLICIT NONE
      REAL(REALK),   INTENT(IN) :: t
      REAL(REALK),   INTENT(IN) :: lambda
      INTEGER(INTK), INTENT(IN) :: v  ! this is the V-point
      INTEGER(INTK), INTENT(IN) :: u  ! this is the rho-point
      REAL(REALK) :: get_F3_LIP3

      INTEGER(INTK) :: k
      REAL(REALK)   :: F1, F2

      k = v-u

      IF (ABS(t) < T_VALUE_THR) THEN

      get_F3_LIP3 = (0.177147D6 / 0.3520D4 * lambda ** 9 - 0.9D1 / 0.10D2   &
        * lambda ** 9 * k ** 7 - 0.729D3 / 0.4D1 * lambda ** 9   &
        * k ** 3 - 0.8748D4 / 0.55D2 * lambda ** 9 * k + 0.35721D5   &
        / 0.160D3 * lambda ** 9 * k ** 2 + 0.3D1 / 0.64D2 * lambda ** 9   &
        * k ** 8 - 0.162D3 / 0.5D1 * lambda ** 9 * k ** 5 + 0.6075D4   &
        / 0.64D2 * lambda ** 9 * k ** 4 + 0.567D3 / 0.80D2 * lambda ** 9   &
        * k ** 6) * t ** 8 + (-0.3D1 / 0.16D2 * lambda ** 7 * k ** 6   &
        - 0.18225D5 / 0.224D3 * lambda ** 7 * k ** 2 - 0.5103D4 / 0.160D3   &
        * lambda ** 7 - 0.243D3 / 0.16D2 * lambda ** 7 * k ** 4 + 0.27D2   &
        / 0.10D2 * lambda ** 7 * k ** 5 + 0.324D3 / 0.7D1 * lambda ** 7   &
        * k ** 3 + 0.2187D4 / 0.28D2 * lambda ** 7 * k) * t ** 6   &
        + (-0.27D2 / 0.5D1 * lambda ** 5 * k ** 3 - 0.972D3 / 0.35D2   &
        * lambda ** 5 * k + 0.9D1 / 0.16D2 * lambda ** 5 * k ** 4   &
        + 0.3645D4 / 0.224D3 * lambda ** 5 + 0.729D3 / 0.40D2   &
        * lambda ** 5 * k ** 2) * t ** 4 + (0.27D2 / 0.5D1 * lambda ** 3   &
        * k - 0.243D3 / 0.40D2 * lambda ** 3 - 0.9D1 / 0.8D1   &
        * lambda ** 3 * k ** 2) * t ** 2 + 0.9D1 / 0.8D1 * lambda

      ELSE

      F1 = (-0.3D1 * exp(-t ** 2 * k ** 2 * lambda ** 2) * lambda ** 2   &
         + 0.4D1 * exp(-t ** 2 * k ** 2 * lambda ** 2) * k * lambda ** 2   &
         - exp(-t ** 2 * lambda ** 2 * (k - 0.3D1) ** 2) * k * lambda ** 2   &
         + exp(-t ** 2 * lambda ** 2 * (k - 0.3D1) ** 2) * k ** 2   &
         * lambda ** 2 - exp(-t ** 2 * k ** 2 * lambda ** 2) * k ** 2   &
         * lambda ** 2) / lambda ** 3 / t ** 2 / 0.4D1 + (-exp(-t ** 2   &
         * k ** 2 * lambda ** 2) + exp(-t ** 2 * lambda ** 2   &
         * (k - 0.3D1) ** 2)) / lambda ** 3 / t ** 4 / 0.4D1

      F2 = (((-k ** 2 + k ** 3 / 0.4D1 + 0.3D1 / 0.4D1 * k)   &
         * erf(-dble(3 * t * lambda) + dble(t) * k * dble(lambda))   &
         + (k ** 2 - k ** 3 / 0.4D1 - 0.3D1 / 0.4D1 * k) * erf(dble(t)   &
         * k * dble(lambda))) / dble(t) + ((-0.1D1 / dble(lambda ** 2)   &
         / 0.2D1 + 0.3D1 / 0.8D1 / dble(lambda ** 2) * k) * erf(-dble(3   &
         * t * lambda) + dble(t) * k * dble(lambda)) + (0.1D1   &
         / dble(lambda ** 2) / 0.2D1 - 0.3D1 / 0.8D1 / dble(lambda ** 2)   &
         * k) * erf(dble(t) * k * dble(lambda))) / dble(t ** 3))   &
         * sqrt(0.3141592654D1)

      get_F3_LIP3 = F1 + F2

      END IF

   END FUNCTION get_F3_LIP3

!-------------------------------------------------------------------------------

   FUNCTION get_F4_LIP3(t,lambda,v,u)

      IMPLICIT NONE
      REAL(REALK),   INTENT(IN) :: t
      REAL(REALK),   INTENT(IN) :: lambda
      INTEGER(INTK), INTENT(IN) :: v  ! this is the V-point
      INTEGER(INTK), INTENT(IN) :: u  ! this is the rho-point
      REAL(REALK) :: get_F4_LIP3

      INTEGER(INTK) :: k
      REAL(REALK)   :: F1, F2

      k = v-u

      IF (ABS(t) < T_VALUE_THR) THEN

      get_F4_LIP3 = (0.3915D4 / 0.64D2 * lambda ** 9 * k ** 4 + 0.190269D6   &
        / 0.3520D4 * lambda ** 9 + 0.15309D5 / 0.80D2 * lambda ** 9   &
        * k ** 2 - 0.33777D5 / 0.220D3 * lambda ** 9 * k + 0.63D2   &
        / 0.20D2 * lambda ** 9 * k ** 6 - 0.2187D4 / 0.16D2 * lambda ** 9   &
        * k ** 3 - 0.351D3 / 0.20D2 * lambda ** 9 * k ** 5 - 0.13D2   &
        / 0.40D2 * lambda ** 9 * k ** 7 + lambda ** 9 * k ** 8 / 0.64D2)   &
        * t ** 8 + (-0.11745D5 / 0.224D3 * lambda ** 7 * k ** 2 - 0.27D2   &
        / 0.4D1 * lambda ** 7 * k ** 4 + 0.351D3 / 0.14D2 * lambda ** 7   &
        * k ** 3 - lambda ** 7 * k ** 6 / 0.16D2 + 0.6561D4 / 0.112D3   &
        * lambda ** 7 * k + 0.39D2 / 0.40D2 * lambda ** 7 * k ** 5   &
        - 0.2187D4 / 0.80D2 * lambda ** 7) * t ** 6 + (0.81D2 / 0.10D2   &
        * lambda ** 5 * k ** 2 + 0.2349D4 / 0.224D3 * lambda ** 5   &
        - 0.1053D4 / 0.70D2 * lambda ** 5 * k - 0.39D2 / 0.20D2   &
        * lambda ** 5 * k ** 3 + 0.3D1 / 0.16D2 * lambda ** 5 * k ** 4)   &
        * t ** 4 + (-0.27D2 / 0.10D2 * lambda ** 3 + 0.39D2 / 0.20D2   &
        * lambda ** 3 * k - 0.3D1 / 0.8D1 * lambda ** 3 * k ** 2)   &
        * t ** 2 + 0.3D1 / 0.8D1 * lambda

      ELSE

      F1 = (0.2D1 * exp(-t ** 2 * k ** 2 * lambda ** 2) * lambda ** 2   &
         + exp(-t ** 2 * k ** 2 * lambda ** 2) * k ** 2 * lambda ** 2   &
         - 0.3D1 * exp(-t ** 2 * k ** 2 * lambda ** 2) * k * lambda ** 2   &
         - exp(-t ** 2 * lambda ** 2 * (k - 0.3D1) ** 2) * k ** 2   &
         * lambda ** 2 - 0.2D1 * exp(-t ** 2 * lambda ** 2   &
         * (k - 0.3D1) ** 2) * lambda ** 2) / lambda ** 3 / t ** 2   &
         / 0.12D2 + (exp(-t ** 2 * k ** 2 * lambda ** 2) - exp(-t ** 2   &
         * lambda ** 2 * (k - 0.3D1) ** 2)) / lambda ** 3 / t ** 4 / 0.12D2

      F2 = (((-k ** 3 / 0.12D2 - k / 0.6D1 + k ** 2 / 0.4D1)   &
         * erf(-dble(3 * t * lambda) + dble(t) * k * dble(lambda))   &
         + (k / 0.6D1 - k ** 2 / 0.4D1 + k ** 3 / 0.12D2) * erf(dble(t)   &
         * k * dble(lambda))) / dble(t) + ((-0.1D1 / dble(lambda ** 2)   &
         * k / 0.8D1 + 0.1D1 / dble(lambda ** 2) / 0.8D1) * erf(-dble(3   &
         * t * lambda) + dble(t) * k * dble(lambda)) + (-0.1D1   &
         / dble(lambda ** 2) / 0.8D1 + 0.1D1 / dble(lambda ** 2)   &
         * k / 0.8D1) * erf(dble(t) * k * dble(lambda))) / dble(t ** 3))   &
         * sqrt(0.3141592654D1)

      get_F4_LIP3 = F1 + F2

      END IF

   END FUNCTION get_F4_LIP3

!-------------------------------------------------------------------------------

END MODULE fefmm_LIP3_Fij
